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ABSTRACT 

Recent observations indicate that star formation occurs only in the molecular phase of a galaxy's 
interstellar medium. A realistic treatment of star formation in simulations and analytic models of 
galaxies therefore requires that one determine where the transition from the atomic to molecular gas 
occurs. In this paper we compare two methods for making this determination in cosmological simu- 
lations where the internal structures of molecular clouds are unresolved: a complex time-dependent 
chemistry network coupled to a radiative transfer calculation of the dissociating ultraviolet (UV) radi- 
ation field, and a simple time-independent analytic approximation. We show that these two methods 
produce excellent agreement at all metallicities > 10 -2 of the Milky Way value across a very wide 
range of UV fields. At lower metallicities the agreement is worse, likely because time-dependent effects 
become important; however, there are no observational calibrations of molecular gas content at such 
low metallicities, so it is unclear if either method is accurate. The comparison suggests that, in many 
but not all applications, the analytic approximation provides a viable and nearly cost-free alternative 
to full time-dependent chemistry and radiative transfer. 

Subject headings: Cosmology: theory — Galaxies: evolution — Galaxies: ISM — ISM: molecules — 
Methods: numerical — Stars: formation 



1. INTRODUCTION 



Resolved observations of nearby galaxies in the last 
decade have unequivocally established that the star for- 
mation rate in a galaxy correlates much more strongly 
with its molecular conten t than with its total gas con- 



tent (W ong fc Bl itz 2002 ; jKennicutt et aL||2QQ7| |Bigiel 



et al." 
Driven 



T pJ8| [Leroy et al.| |:M)8| |BIanc et affg )09 ) . 

m by these observations, analytic m odels (Blitz & 



Rosolowsky 2006; Krumh olz et al.Jl2009c|), semi-analytic 



models (Obreschkow fc Rawlings 2009; Obreschkow et al. 

2010j |Dutton et al.||2010|), and Simula- 



2009 



Fu et al, 



tions (FeTu pessy et al. 2006; 
Pelupes sy fc PapadopouT* 



Gnedin fc Kravtsov| 2010a|o'| ) have alf beg 
porate models tor the H I to H2 transitio 
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2008 
2009 



;un to mcor- 
2 transition into their 
recipes for star formation. The inclusion of molecule for- 
mation has allowed these models to solve a number of 
problems that appeared in earlier models that omitted 
molecules. For example in the high redshift universe, 
absorption line surveys indicate that the star formation 
rate in damped Lyman a systems must be significantly 
below what one would expect given their H I column 
densities if the star formation rate were correlate d with 
total gas surface density as it is i n nearby galaxies (IWolfe 



& Chen 2006 



Wild et al. 2007). This can be explained 



by the tact that these systems are usually quite metal 
poor, which causes them to be systematically deficient 
in molecular hydrogen and thus in star formation com- 
pared to local metal-rich galaxies of similar total gas col- 
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umn d ensity ( Krumholz et 



2010b). Similarly, in the loca. 



al.|2009a|[Gnedm 
>cal universe, moc 



& Kravtsov 
dels incorpo- 
rating molecule formation are able to explain why some 
galaxies that undergo H I stripping in a cluster environ- 
ment have their star formation quenched, while others 
do not. Those that lose their outer gas disks but re- 
tain high column density, molecule-rich centers continue 
to form stars, while those that lo se their inner gas disk s 
cease to form molecules and stars ( Fumagalli et al.|2 009). 

While many of these authors have compared their 
models or simulations to observed molecular fractions in 
nearby galaxies, these observations cover only a limited 
dynamic range in metallicity, radiation environment, and 
other quantities that may be relevant to the H I to H2 
transition. We cannot yet test the models in more ex- 
treme environments with observations, but we can check 
for consistency between the models. Such checks provide 
at least a minimum level of confidence in extrapolating 
the model predictions beyond the range of environments 
in which they have been calibrated. Cross-comparisons 
are particularly important given the range of complex- 
ity of the models, which go from steady-state analytic 
approximation formulae to sophisticated time-dependent 
radiation plus chemistry modules. Only the former are 
suitable for implementation in semi-analytic models, and 
thus such a test is necessary to ensure that semi-analytic 
models that rely on them are able to faithfully reproduce 
the results of numerical simulations. 

Our goal in this paper is to provide such an inter- 
model com£aris^rL__We_Qse a sui te of numerical mod- 
elsby I Gnedin fc Kravtsov] (|2010b| hereafter GK10; see 



2.1 



lor a detailed description), and we compare them 
to on e of the most commonly-used analytic approxima- 
tions ([Krumh olz et al.||2008[ |2009bj |McKee fc Krumholz 



2010 see \T71 tor a detailed description) tor the H I to 



H2 transition. Both of the analytic and numerical models 
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have been tested against observations of local galaxies, 
and have been found to provide a good description of the 
molecular content in themj^] This comparison will allow 
us to establish under what circumstances the numerical 
and analytic models agree and thus can be applied with 
reasonable confidence, and to gain physical insight into 
the H2 formation process by examining cases where they 
differ. An additional practical benefit is that, in cases 
where the analytic models are able to do a reasonable 
job reproducing the numerical ones, they may provide an 
approximation that can be used in simulations at greatly 
reduced computational cost. 

The structure of the remainder of this paper is as fol- 
lows. In § [2] we summarize the numerical and analytic 
models for the H I to H2 transition. In § [3] we compare 
the models in idealized simulations in which the radia- 
tion field and metal content are held fixed, while in 
we consider simulations in which the radiation field an 
metal content are determined in a fully self-consistent 
manner. Finally in § [5] we discuss and summarize our 
conclusions. 

2. SUMMARY OF THE METHODS 

Here we summarize how we estimate the H2 mass 
fraction in every com putational cell of our test simu- 
lations. Section |2.1| describes the ti me-de pendent nu- 
merical chemistry method, and Section [Z2| describes the 
analytic model and how we apply it to the simulations. 

2.1. The Numerical Model 

Simulations used in this paper are fully described in 
GK10. Here we only remind a reader that a set of "fixed 
ISM" simulations was perform ed with the Adaptive Re- 
finement Tree (ART) code ([K ravtsov 1999 Kravtsov 



et 



aL] |2002| |Rudd et aL]|20()8p . Each simulation was 
rted from a z = 4 output of a fully self-consistent cos- 



started from a z = 4 output of a fully self-consistent cos- 
mological simulation of a Milky Way progenitor galaxy. 
However, in each fixed ISM simulation the dust-to-gas ra- 
tio and the interstellar radiation field were fixed to pre- 
defined, constant in space and time values, effectively 
turning a cosmological simulation into a highly realistic 
simulation of an isolated galaxy with a live dark matter 
halo and natural accretion history. 

Each fixed ISM simulation was continued for 600 Myr, 
long enough to fully establish a new ISM structure and 
eliminate the memory of the initial conditions. The dark 
matter mass resolution in the simulation is 1.3 x 10 6 M 
and the spatial resolution of the highest resolved region of 
the simulation (that includes the whole of the simulated 
galaxy disk) is 65 pc. That translates into the mass res- 
olution in the gas of between ~ 10 3 M and ~ 10 6 M , 
depending on the value of the gas density. 

All fixed ISM simulations include gas dynamics, 
star formation, and the time-dependent and spatially- 
inhomogeneous 3D radiative transfer of ultraviolet and 
ionizing radiation from individual stellar particles that 
formed during the course of each simulation. The simu- 
lations incorporate a non-equilibrium chemical network 
of hydrogen and helium and non-equilibrium cooling and 

5 In contrast, we do not compare to the H^midplane pressure 
correlation model of |Blitz Rosolowsky (2006), because recent 
observations indicate that it breaks down on the sub-kpc scales 
resolved in the simulations (Fumag alli et al.|2010| ). 



heating rates, which make use of the local abundance of 
atomic, molecular, and ionic species and UV intensity. 
This network includes formation of molecular hydrogen 
both in the primordial phase and on dust grains. The 
abundances of the relevant atomic and molecular species 
are therefore followed self-consistently during the course 
of the simulation. The heating and cooling terms in the 
equation for the internal energy include all of the terms 
normally included in the simulations of first stars and in 
the ISM models, including cooling on metals. A complete 
set of chemical reactions and the values for the adopted 
reaction and cooling/heating rates are presented in the 
appendix of GK10. 

In addition to the simulations with the fixed ISM prop- 
erties, we also use a fully self-consistent cosmological sim- 
ulation of several galaxies in a cosmological volume with 
25h~ 1 Mpc on a side. These simulations will be fully de- 
scribed elsewhere (M. Zemp et al, in preparation), but in 
physical modeling they are identical to the simulation of 



|Gnedin fc Kravtsov] pO 10a ); the only difference from the 
early simulation is a larger box size and slightly different 
values of cos mological parameter s, based on the WMAP5 
cosmology ( |Dunkley et al.|[2008| ). The spatial and mass 
resolution of the cosmological simulation is matched ex- 
actly to the spatial and mass resolution of the simulation 
with the fixed ISM. Due to computational expense, the 
fully self-consistent cosmological simulation was not con- 
tinued beyond z = 3. For comparison with the analytical 
model, we use several of the model galaxies from that 
simulation that differ widely in their ISM properties. In 
model galaxies from a full cosmological simulation the 
properties of the ISM (dust-to-gas ratio, gas metallicity, 
the interstellar radiaton field, the gas density distribu- 
tion, etc.) vary inside the galaxy in a manner that is 
controlled by their prior cosmic history and the present 
cosmological environment. 

Finally, we call the reader's attention to two details of 
the microphysical model of H2 formation in the simula- 
tions that will be important for our comparison. First, 
since the H2 destruction rate depends on the amount of 
shielding, one must estimate a column density for each 
cell. In the simulations one does this using the Sobolev 
approximation: in each cell one computes the local den- 
sity scale height h = p/|Vp|, and then takes the column 
density to be U = ph. Comparison of this approxima- 
tion with ray-tracing done in post-process shows that it is 
generally quite accurate. Second, one must increase the 
rate coefficient for H2 formation by a clumping factor to 
account for density inhomogeneities that are unresolved 
on the computational grid. Since the H2 formation rate 
varies as the square of density, these inhomogeneities in- 
crease the overall rate. The simulations use a clump- 
ing factor of 30, which is chosen to give a good match 
between the simulations and FUSE observations of the 
Milky Way and Magellanic clouds. 

2.2. The Analytic Model 

The analytic model to whic h we compare th e simu - 
lations is based on the wor k of Krumholz et al.| ([2008), 
Krumholz etaT] ( |2009b| ) , an dllVlcKee & Krumholz] pOTOp , 
hereafter known as the KMT model. KMT consider an 
idealized spherical cloud immersed in a uniform, isotropic 
Lyman- Werner band radiation field. They then solve the 
coupled problems of radiative transfer and H2 formation- 



3 



dissociation balance under the assumption that the cloud 
is in steady state. They do not co nsider gas-phase reac- 
tions (though see the Appendix of |McKee fc Krumholz 
2Q1Q| for a discussion of how these might be included). 
The solution is 



/h 2 



1 + 0.25s' 



where /h 2 is the H 2 mass fraction in the cloud, 
ln(l + 0.6x + 0.01x 2 ) 



0.6r c 



X = 71 



Cfd-21 

n. 



16.5 



(1) 

(2) 
(3) 



r c is the dust optical depth of the cloud, (Jd -21 is the 
dust cross-section per H nucleus to 1000 A radiation, 
normalized to a value of 10 -21 cm -2 , 7£_i6.5 is the rate 
coefficient for H2 formation on dust grains, normalized to 
the M ilky Way value of 10" 16 ' 5 cm 3 s" 1 (IWomre et al. 
|2008[ ) , Gq is the ambient UV radiation field intensity, 
normalized to the |Draine (|1978[ ) value for the Milky Way, 
and nn,o is the volume density of H nuclei in units of 
cm -3 . Note that cr d and 1Z are both proportional to the 
dust to gas ratio, so their ratio is independent of the 
metallicity. 

In order to apply this analytic model to the simula- 
tions, we must have a way of computing the two di- 
mensionless numbers \ an d t c on which fn 2 depends. 
Our estimate of r c is straightforward. Since the reso- 
lution of the simulations is well-matched to individual 
molecular clouds, we estimate the local column density 
E from the Sobolev ap prox imation exactly as the simu- 
lations do (see Section 2.1). The dust optical depth is 
simply r c = Ead//iH, where ad is the dust cross section 
per H nucleus and /i H = 2.3 x 10~ 24 g is the mean mass 
per H nucleus. We take the cross section to have a value 
&d,-2i = Z' ', where Z' is the metallicity normalized to 
the Milky Way value. For Z' we use the local metallicity 
in the cell for which we are computing r c . 

The scaled radiation field x requires a bit more care. 
In Section |3] we discuss simulations where the UV field 
is fixed to a specified value rather than determined self- 
consistently. In that case we compute x directly from 
equation O) using the UV field, gas volume density, dust 
cross section, and rate coefficient adopted in the simula- 
tion. The dust cross section is again &d-2i = Z while 
the rate coefficient is 1Z-iq^ = 30Z' . This value of 1Z 
includes the factor of 30 enhancement due to unr esolv ed 
clumping adopted in the simulations (see Section |2.1[ ). 

In contrast, in Section [4] we compare to simulations 
where the UV field is determined self-consistently and 
varies from point to point. In this case we can proceed 
in two ways. First, we can use the UV field computed in 
the simulations in every cell, in analogy to the first case. 
Second, and more interestingly, we can omit this infor- 
mation. KMT show that, if the interstellar medium is 
in a self-consistently determined two phase equilibrium, 
the ratio Go/nn,o should take on a characteristic value 
such that 

1 + 3.1Z /0365N 



where Z' is the metallicity relative to solar units. This 
is not likely to hold cell-by-cell at every time step, but 
KMT argue that it should hold on average. The ad- 
vantage of this approach is that using this value of x 
can we can compute the H2 fraction using only the local 
metallicity and the column density as determined by the 
Sobolev approximation. Thus this method can be used in 
simulations or semi-analytic models that do not explic- 
itly include either chemistry or radiative transfer. Since 
the radiative transfer and chemistry are often the most 
computationally costly parts of a simulation, it would be 
extremely useful to be able to roughly approximate the 
results of a full chemistry and radiation simulation with 
such a simple analytic model. 

3. SIMULATIONS WITH FIXED ISM PROPERTIES 

In order to gain a better understanding of how the 
model and simulations compare, we begin by consid- 
ering a set of simulations in which the radiation field 
and metallicity are specified by hand to have a particu- 
lar value, rather than being determined self-consistently. 
The grid of metallicities runs from Z' = 10~ 3 - 3, and 
the grid of UV fields runs from G' = 0.1 — 100. From 
each simulation we extract all computational cells that 
are located in the disk of the simulated galaxy. For each 
computational cell we record its total hydrogen density, 
temperature, metallicity, gas-to-dust ratio, local value of 
the interstellar radiation field, and abundances of ion- 
ized, atomic, and molecular hydrogen. 

To quantify the level of agreement between the simu- 
lations and the analytic model, in every computati onal 
cell we use the analytic model described in Section |2.2 
to compute a predicted H2 fraction /H 2 ,anaiyt- We com- 
pute the absolute difference between this and the value 
obtained in the simulation /H 2 ,sim? 



A/h 2 = |/H 2 ,analyt ~ /h 2 ,s 



(5) 



X»3.1 



4.1 



(4) 



because the star formation rate depends on the absolute 
value of the molecular abundance, not on relative one. 
For example, regions with /h 2 <C 1 may show a large 
relative error between the two models, but those regions 
would not contribute to star formation rate in a galaxy 
or on a particular spatial scale. 

Figure [l] shows /(< A/h 2 ), the fraction of the mass in 
the extracted galaxy disk cells (not in the entire cosmo- 
logical simulation) for which the difference is less than 
A/h 2 . In this plot perfect agreement would appear as 
a function of constant value unity, while complete dis- 
agreement would be a function of constant value zero. 
The Figure indicates that the agreement is extremely 
good for all UV field strengths at metallicities Z' > 10 -2 . 
In all of these cases the disagreement is essentially zero 
for at least half of the mass, and is less than 25% for 
> 75% of the mass. The Milky Way-like case Z' = 1 
and Gq = 1 shows near-perfect agreement between the 
simulations and the analytical model. In contrast, the 
agreement is quite poor for the lower metallicity cases, 
particularly for the combination of low metallicity and 
weak radiation field. 

In order to understand why some simulations agree and 
others do not, we construct bins of /H 2 ,sim running from 
- 0.2, 0.2 - 0.4, 0.4 - 0.6, 0.6 - 0.8, and 0.8 - 1.0, and 
similarly for /H 2j anaiyt- We then compute what fraction 
of the total mass in the simulation falls into a given bin 
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Fig. 1. — Cumulative fraction of the total simulation mass /(< A/h 2 ) f° r which the difference between the simulation and analytically- 
predicted values of /h 2 is less than A/h 2 - Columns show simulations of fixed metallicity relative to Milky Way, Z x , as indicated at the 
top of each column, and rows show simulations of fixed UV radiation field relative to Milky Way, Gq, as indicated at the left of each row. 
Dashed horizontal lines show a value of 1.0, indicating perfect agreement. An X in a box indicates that no simulation for that value of Z' 
and Gq is available. 



0f/H 2 . sim versus /H 2 ,anaiyt- If agreement were perfect all 
the mass would fall into the same bin on either axis, i.e. 
all the mass with /H 2 ,sim in the range — 0.2 would also 
have /H 2 ,anaiyt m this range, and so forth. 

We plot the actual distribution of mass in bins in Fig- 
ures [2] and [3J The former figure shows a linear color 
scale, indicating where the bulk of the mass lies, while 
the latter shows a logarithmic scale, intended to empha- 
size areas of disagreement even if they contain relatively 
little mass. The figures reveal some interesting system- 
atic differences between the simulations where the an- 
alytic model agrees well with the numerical results, at 
Z' > 10 -2 , and the ones where it does not, at Z' < 10 -2 . 
First, in the higher metallicity cases the majority of the 
mass in both the simulations and the analytic models 
tends to be mostly molecular or mostly atomic; bins of 
intermediate molecular fraction are sparsely populated. 
In contrast, at lower metallicity a significant fraction of 
the mass is in the intermediate molecular fraction regime, 
and here the analytic model tends to systematically over- 
estimate the molecular fraction compared to the simula- 
tion. 

Although the cause of the discrepancy is not entirely 
clear, we regard the most likely explanation to be that 
the steady-state assumption in the model is breaking 
down at metallicities Z' < 10 -2 . Consider a compu- 
tational cell of volume density nn that is predominantly 



atomic. If the dissociation rate is negligible, the time 
required for that cell to convert to molecular is the ra- 
tio of the H I volume density to the volumetric rate of 
conversion of H I to H2 

* H ° = 2^R = lln ^ Z " 1 Myr ' ( 6) 

where we adopt the rate coefficient 7£_i6.5 = 30Z' used 
in the simulations. 

Figure [4] illustrates the implications of this timescale 
by comparing one simulation where the analytic model 
agrees very well with the simulation H2 fraction [Z 1 = 1, 
Gq = 1) to another where they agree poorly (Z 1 = 0.003, 
Gq = 1). In the Milky Way metallicity case the locus 
where the analytic model predicts 50% H2 fraction inter- 
sects the region of parameter space occupied by the sim- 
ulation points near nn ~ 10 cm -3 , corresponding to an 
H2 formation time tn 2 ~ 2 Myr at Z' = 1. This is short 
compared to any galactic timescale, and thus the instan- 
taneous equilibration approximation is an excellent one. 
On the other hand, in the low metallicity simulation the 
locus of intersection is at nn ~ 100 cm -3 , corresponding 
to tn 2 ~ 100 Myr at Z' = 0.003. In this case a parcel 
of gas may be compressed in a spiral arm to the point 
where one might expect it to become H2-dominated were 
it able to reach steady state, but before it has a chance to 
do so it leaves the spiral arm and its density and column 
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0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 

f|H2,sim ^H2,sim ^H2,sim ^H2,sim ^H2,sim W.sim ^H2,sim ^H2,sim 

Fig. 2. — In each panel the colors pixels indicate the fraction of the mass of that simulation in the given bin of /H 2 ,sim versus /H 2 ,analyt 5 
divided by the fraction of mass in the most massive bin for that simulation, so that the bin in each panel that contains the most mass is 
always black. See text for details. Perfect agreement would correspond to all pixels lying along the dashed one-to-one lines through the 
center of each panel. Simulations are labeled with values of Z' and G' Q as in Figure^ 



Z 1 = 0.001 Z 1 = 0.003 Z 1 = 0.01 Z 1 = 0.03 Z 1 = 0.1 Z 1 = 0.3 Z 1 = 1.0 T = 3.0 




0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.0 0.5 



1 H2,sim 'H2,sim 1 H2,sim 1 H2,sim 1 H2,sim 1 H2,sim 1 H2,sim 1 H2,sim 

Fig. 3. — Same as Figure [2] but using a logarithmic rather than a linear scale to emphasize areas of disagreement. 
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Fig. 4. — Each point shows the value of nn and r c in a compu- 
tational cell in the simulation with the indicated values of Z' and 
Gq. (To minimize file size we show only a subset of cells.) Points 
are color-coded by molecular fraction in the simulation, with black 
indicating /H 2 ,sim < 0.1, red indicating 0.1 < /H 2 ,sim < 0.5, green 
indicating 0.5 < /H 2 ,sim < 0.9, and blue indicating 0.9 < /H 2 ,sim- 
In contrast, the solid line separates points with /H 2 ,analyt < 0.5 
(below the line) from those with /H 2 ,analyt > 0.5 (above the line). 

density fall. The analytic model identifies such regions 
as predominantly molecular because it ignores the time 
dependence. 

Before concluding that the simulations are correct and 
the model erroneous in this case, however, we must add 
a final caution. The timescale we compute is inversely 
proportional to the subgrid clumping factor assumed in 
the simulations. While this number has been estimated 
for metallicities down to Z' « 0.2 in the Small Magel- 
lanic Cloud, we have no direct knowledge of its value in 
lower metallicity galaxies. Nonetheless, we tentatively 
conclude that the time-dependence of the H2 fraction is 
not a major effect at metallicities Z' > 10 -2 , so a time- 
independent model is adequate in this regime. 



4. 



SIMULATIONS WITH SELF-CONSISTENTLY 
DETERMINED RADIATION FIELDS 



We now proceed to compare the KMT model to simula- 
tions that include an explicit, self-consistent calculation 
of the local radiation field, and use this local radiation 
field to compute the H2 dissociation rate. For this com- 
parison we s elect three galaxies from the simulations de- 
scribed in § |2.1| at redshift z — 3, chosen to represent 
a wide range of galactic types and environments. Table 
[l] lists some of the basic properties of the three galax- 
ies. Galaxy 1 is a typical z — 3 star- forming galaxy, 
with a high radiation field, relatively high metallicity, 
and generally high gas densities; by z = it is expected 
to evolve into a Milky Way type galaxy. In contrast, 
galaxy 2 has just started its first starburst. As a result, 
it has a fairly strong radiation field, but it has not yet 
produced many metals and thus has a very low metallic- 
ity. Finally, galaxy 3 is a dwarf, SMC-like galaxy that 
has been forming stars at a low rate, and remains mostly 
gaseous. As a result is has the lowest radiation field of 
any of the three, but a fairly high metallicity. 

Figure [5] shows, for each of the three sample galax- 
ies, the cumulative fraction of the simulation mass for 
which the difference between the simulation and the 
KMT model is less than A/h 2 - The dashed line cor- 
responds to the case where we compute A/h 2 using our 
knowledge of the local radiation field, while the solid line 
corresponds to adopting the mean radiation field given 
by equation Q. The latter is the model that can be ap- 
plied in a simulation that does not include a radiation 
and chemistry module. 

The first thing to notice is that, for all three galaxies, 
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Fig. 5. — Same as Figure ^ but for the three sample galaxies 
with self-consistently computed radiation fields and H2 fractions. 
Dashed lines show the result when we compute /h 2 m the KMT 
model using the local radiation field in each cell, while solid lines 
show the result when we compute /h 2 in the KMT model with- 
out using any information about the local radiation field, and we 
instead use the mean radiation field given by equation E|. 

the KMT model using the mean radiation field (solid 
lines) agrees very well with the simulation result. In all 
cases the difference in H2 fraction between the model and 
simulation is less than 10% for ~ 80% of the mass, and 
is less than 20% for ~ 90% of the mass. This indicates 
that we are able to reproduce the simulation results very 
well, across a wide range of galaxies, using a model that 
requires no knowledge of the local radiation field. This 
is highly encouraging for numerical simulations and ana- 
lytic and semi-analytic models, since it suggests that the 
results of a full radiative transfer and chemistry solve 
can be approximated reasonably well at far less compu- 
tational cost. 

Surprisingly, the agreement between the simulations 
and the KMT model actually worsens if we do use our 
knowledge of the local radiation field. To help under- 
stand why this is, in Figures [6] and [7] we show plots anal- 
ogous to those in Figures [2| and |3| for the self-consistent 
case. As in Figures [2] and |3f we take each simulation and 
divide the mass into bins of /H 2jS im and /H 2 ,anaiyt> and 
plot how much mass falls into each bin; perfect agreement 
would then appear as all colored pixels falling along the 
diagonal 1 — 1 line. The plots show that, in the case 
where we adopt the mean radiation field of equation Q, 
most of the mass is indeed in pixels along the 1 — 1 line, 
and mass that is off this line is scattered fairly symmet- 
rically about it. In contrast, if we use the instantaneous, 
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TABLE 1 

Summary of Galaxy Properties 



Galaxy 


Halo mass (Mq) 


Stellar mass (Mq) 


Gas mass (Mq) 


(/hi) 


</h 2 > 


(Z f ) 


( g o)m 


(G' )v 


1 


4.8 x 10 11 


4.4 x 10 10 


4.6 x 10 10 


0.09 


0.02 


0.5 


1,100 


92 


2 


2.3 x 10 10 


5.3 x 10 7 


5.9 x 10 9 


0.37 


0.09 


0.01 


160 


17 


3 


2.9 x 10 10 


1.0 x 10 9 


4.5 x 10 9 


0.26 


0.15 


0.18 


39 


7.4 



Note. — Col. 5-8: Mass- weighted average H I fraction and H2 fraction within the virial radius of a galaxy, 
relative metallicity of neutral gas, and relative radiation intensity in the neutral gas (mass and volume weighted). 
Note that (/hi) + (/h 2 ) < 1 because the gas mass budget includes ionized hydrogen as well. 




0, 

0.0 I 

0.8 - 
> 0.6 - 
0.4 '- 
0.2 , 



0.0 

0.0 0.2 0.4 0.6 0.8 



0.8 
f 0.6 
J 0.4 
0.2 
0.0 

0.8 

f °- 6 
*5 0.4 
0.2 



- 0.6 
0.4 
0.2 



1 



0.0 0.2 0.4 0.6 O.f 



1 1 

0.0 0.2 0.4 0.6 0.1 



0.0 0.2 0.4 0.6 O.f 



Fig. 6. — Same as Figure [2] but for the three sample galaxies 
with self-consistently computed radiation fields and H2 fractions. 
The left column shows the results using the mean radiation field 
given by equation El, while the right column gives the results using 
the local radiation neld in each cell computed in the simulation. 

known radiation field value in each cell, we find that the 
KMT model systematically overestimates the molecular 
fraction relative to the simulation. This effect is small- 
est in galaxy 1, which has the highest metallicity, and is 
larger in galaxies 2 and 3, which have lower metallicity. 

This is very likely another version of the timescale ef- 
fect we identified in Section 03 The timescale to reach 
equilibrium in the lower metallicity simulations is not 
much smaller than the timescale over which radiation 
fields and densities fluctuate, and thus the relevant value 
of x is a time- averaged rather than an instantaneous one. 
The instantaneous value of G' is usually smaller than its 
time-average (because there are few regions of intense 
radiation and many regions of little radiation), and so 
using this instantaneous value leads to an overestimate 
of the H2 fraction in those cases where the metallicity is 
low enough to make the equilibration time long. It is in- 
teresting, however, that this effect seems to be mitigated 
if we simply adopt the mean characteristic radiation field 
expected for two-phase equilibrium. In this case we ob- 
tain good agreement even with for the lowest metallicity 
self-consistent galaxy, which has (Z f ) = 0.01. 

5. SUMMARY 



Fig. 7. — Same as Figure [6] but using a logarithmic rather than 
a linear scale to emphasize areas of disagreement. 

In this work we compare two approaches to the prob- 
lem of determining where the ISM in galaxies transitions 
from H I to H2. One is a highly- accurate time-dependent 
chemistry and radiative transfer method implemented 
in a large-scale cosmological simulation code, while the 
other is a simple analytic approximation suitable for im- 
plementation in semi-analytic models or simulation codes 
that to not include chemistry or radiative transfer. Our 
comparison shows that the two models agree extremely 
well whenever the metallicity is of order 1% of the So- 
lar value or more. We find this agreement in both ex- 
perimental galaxies with fixed metallicity and radiation 
field, and in galaxies whose metallicity and radiation field 
are computed in a fully self-consistent manner. The two 
models diverge at even lower metallicities, and we inter- 
pret this as being a result of a failure of the assumption 
of chemical equilibrium in the analytic model. At such 
low metallicities, the time required to convert a given 
fluid element from H I to H2 is comparable to a galactic 
rotation period even in the absence of dissociating radi- 
ation, and the equilibrium model overestimates the H2 
fraction in this case. 

Nonetheless, the excellent agreement we find at metal- 
licities of 1% Solar or more suggests that the analytic 
model has its uses. In simulations, it is capable of repro- 
ducing the results of the much more complex and accu- 
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rate time-dependent chemistry and radiation module at 
essentially zero cost in terms of both computational time 
and programming effort. The analytic model also allows 
analytic or semi-analytic models to incorporate the ef- 
fects of the H I to H2 transition in a way that faithfully 
reproduces the results of simulations. For these reasons 
it is a valuable tool for cosmology. 
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